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A new projection method, implemented in higher-dimensional reciprocal space, is 
developed to compute quasicrystals. The approach enables us to represent quasicrys- 
tals as periodic structures in higher-dimensional space. A proposed projective ma- 
trix can project a higher- dimensional periodic structure into a quasicrystal in physical 
space. Compared with the traditional crystalline approximant method, the projection 
method overcomes the restrictions of the Simultaneous Diophantine Approximation, 
and can use periodic boundary conditions conveniently. Meanwhile, the proposed 
method can efficiently reduce the computational complexity through implementing 
in a unit cell and using pseudospectral method. By applying the projection method 
to the Lifshitz-Petrich model, we can compute quasicrystals rather than crystalline 
approximants, and maintain the rotational symmetry accurately. More significantly, 
the algorithm can calculate the free energy density to high precision. 
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I. INTRODUCTION 



As early as the 1890s, the periodic structures (crystals) in three dimensions were deter- 
mined by 230 space groups based on periodicity, and then the classical crystallography was 
completed, in which the allowed rotational symmetry is only 1-, 2-, 3-, 4-, 6-fold symmetry. 
Both the structure determination and the study of physical properties are based on the 
periodicity which allows the study to be simplified to a unit cell. However, in the 1980s, a 
forbidden 5-fold symmetry electron diffraction pattern was discovered by Shechtman et al.^ 
in a rapid cooled Al-Mn alloy. Later, the term "quasicrystals" appeared for the first time^ 
to describe the non-conventional ordered structures. In an idealized description, quasicrys- 
tals have quasiperiodic, rather than periodic, translational order with non-crystallographic 
symmetry. In the early theory of tilings, the discovery of aperiodic tiling with 5-fold symme- 
try of the plane by Penrose^ showed that such well-ordered systems were mathematically 
possible. Since the original discovery, hundreds of quasicrystals have been reported and con- 
firmed in metallic alloys with 5-, 8-, 10-, 12- fold orientational symmetry^. Shechtman also 
win the Nobel Prize in Chemistry in 2011 for the discovery of quasicrystals. Two decades 
after the first discovery of quasicrystals in metallic alloys, several soft quasicrystals have 
been found in soft matter systems. A number of materials successively joined the family 
of soft quasicrystals, including supramolecular dendrimer or tetraphilic liquid crystals^^, 
ABC star terpolymers^, surfactant-coated metallic nonoparticles^, and colloidal copoly- 
mer micelles in solution 1 -^. Recently the natural quasicrystals have also been discovered in a 
khatyrkite-bearing sample^. To date, quasicrystals are confirmed as a new class of ordered 
structures. 

Two kinds of theoretical approaches have been developed to study quasicrystals, or more 
generally, aperiodic crystals. The first method is an approximate approach which studies 
crystalline approximants rather than quasicrystals. The crystalline approximants are peri- 
odic structures in which the arrangements of lattices closely approximate the local structures 
in quasicrystals. Many crystalline approximants related to quasicrystals have also been dis- 
covered^. These approximants may play an important role in describing the local structures 
of quasicrystals, their formation, stability, and physical properties. The second approach is 
a direct method to study the quasicrystals or aperiodic crystals in the hyperspace, called the 
higher- dimensional approach. In this method, a quasiperiodic structure may be viewed as 
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a periodic structure by extending it into a higher-dimensional space. Its symmetries can be 
expressed in terms of the conventional point groups and space groups of higher-dimensional 
periodic crystals^^. A simple and famous example is the one- dimensional Fibonacci se- 
quence. The sequence is formed by iterative application of the substitution rule : L — > LS, 
S — > L to the two-letter alphabet {L,S}. Then the n-th Fibonacci number F n satisfies 
the recurrence relation F n+ i = F n + The ratio lim^oo F n /F n ~\ = r = (a/5 + l)/2 

is irrational, and is the golden ratio. Therefore the Fibonacci sequence is a quasiperiodic 
series. From another point of view, as Fig. [T] shows, the one-dimensional Fibonacci sequence 
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FIG. 1: Fibonacci sequence obtained as the projection of the square lattice inside the 
window whose width \/2a is shown between broken and solid lines with a slope of 

a = 1/r = (a/5 - l)/2. 

can be obtained as a projection from a two-dimensional square-lattice of points. If the 
slope a is rational, the projected series becomes a periodic chain. In the higher-dimensional 
description, quasiperiodic structures result from irrational physical-space cuts of appropri- 
ate periodic hypercrystal structures. The higher-dimensional approach reveals the hidden 
structural correlations. To implement the higher- dimensional approach in the direct space 
one must know the discrete lattice arrangements of higher- dimensional periodic structure. 
The embedded spaces of (i-dimensional quasiperiodic structures are abstract spaces whose 
dimensions are more than three. The dimensions of the embedded space are dependent on 
the symmetry of the quasicrystal (d > 1) 17 . For example, the quasicrystals with 5-, 8-, 
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10-, and 12-fold symmetry need to be embedded into four- dimensional space. While for 
the quasiperiodic structures with 7-, 9-, 18-fold symmetry, the dimension of the embedding 
spaces increases to six. For other symmetries, embedding spaces with even higher dimension 
will be needed^. Therefore it can make the method difficult to implement. Another dis- 
advantage of applying the higher-dimensional approach in direct space is that in practical 
problems one is required to compute the continuous quasiperiodic distribution rather than 
the discrete quasiperiodic lattice. A convenient way to describe the quasiperiodic structures 
is in the higher-dimensional reciprocal spaced By redefining point-group symmetry and 
space-group symmetry in terms of gauge functions, a broader Fourier-space crystallogra- 
phyESHSl was developed to describe quasicrystals. 

In this paper, we focus on the development of numerical methods for generating qua- 
sicrystals based on the idea of Fourier-space crystallography^^. Quasicrystals are one 
kind of space-filling structures. The traditional numerical method uses a periodic structure 
to approximate the quasicrystal^. In other words, those methods compute a crystalline 
approximant. A natural expectation is that the obtained approximants should approximate 
quasicrystals as the computational box goes to infinity. The advantage of the approximate 
method is that using the periodic boundary conditions is convenient. However, the ap- 
proach cannot obtain quasicrystals exactly unless the computational box is the whole space 
because of the restriction of the Simultaneous Diophantine Approximation (SDA)^. We 
will further explain this point in Sec.|TTJ We then will provide a new numerical method 
based on the higher-dimensional reciprocal space description 119 1 20 1 to calculate quasicrystals 
rather than crystalline approximants. The proposed numerical method is able to overcome 
the restriction of SDA and compute the energy density to high accuracy without periodic 
approximation. 

Although our method is applicable to any model including quasicrystals, for definiteness, 
we use the Lifshitz-Petrich model 122 to demonstrate our algorithm. Before we go further, 
a short introduction to the Lifshitz-Petrich model is necessary. Lifshitz-Petrich model is a 
mean-field theory, and its free energy density functional is 

FW,y)\ = i J dxdy {^[(V 2 + 1)(V 2 + g 2 )0] 2 - |</> 2 - |</> 3 + \<f). (1) 

In Eqn. ([!]), (j)(x, y) is the order parameter. V is the system volume, q is an irrational number 
depending on the symmetry, e is the reduced temperature, c > is an energy penalty to 
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ensure that the principle reciprocal vectors of structures is located on |k| = 1 and |k| = q. 
a > is a phenomenological parameter. For quasicrystals the system volume V should go to 
infinity since quasicrystals are the space-filling structures without periodicity. For ordered 
phases including periodic and quasiperiodic ones, the order parameter cff is the minimum 
of the free energy density functional, which means 

8F 

In order to find the equilibrium state, we introduce a relaxational dynamical to minimize 
the Lifshitz-Petrich energy functional which yields 

f = ~ = -c(V 2 + 1) 2 (V 2 + m + s0 + - <p\ (3) 

We choose a semi-implicit scheme to solve the dynamical equation (|3j): 

^t+At ~ <k) = e<t>t ~ c(V 2 + 1) 2 (V 2 + q 2 ) 2 <Pt + At + «(0 2 ) t - (<P% (4) 

At is the time step size. The specific implementation of the semi-implicit method for different 
numerical methods will be discussed in Sec.HTl 



II. NUMERICAL METHODS 

A. Crystalline Approximant Method (CAM) 

1. Method Description of CAM 

Numerical methods designed for periodic structures have been used to study quasicrystals 
approximately^. Here we call this method the "crystalline approximant method (CAM)". 
In order to describe this method, a brief introduction of numerical methods for treating 
periodic structures is necessary, more details can be found in [21]. For any <i-dimensional 
periodic function /(r), r G 1Z , the repeated structural unit is called a unit cell. A primitive 
unit cell, described by d vectors ei, e 2 , . . . , e^, has the smallest possible volume. The Bravais 
lattice vector is then defined by 

R = he t + l 2 e 2 H h lde d , (5) 
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where 1 = (h, h, ■ ■ ■ , Id) is a d-dimensional vector with components U G Z. For any R in the 

Bravais lattice, the structure is invariant under a lattice translation, i.e., /(r + R) = /(r). 

Given the primitive vectors (e 1; e 2 , . . . , e d ), the primitive reciprocal vectors eJj, . . . , e^) 
satisfy the equation 

e; • e* = 2tt(%. (6) 

The reciprocal lattice vector is then specified by 

k = he{ + k 2 e* 2 + --- + k d e* d , (7) 

where k{ G Z. One of the most important properties of the reciprocal lattices is that plane 
waves {e* k r } form a set of basis functions for any function with periodicity of the lattice. 
The periodic function /(r) on the Bravais lattice can be expanded as 

/(r) = (8) 

k 

For periodic structures, the reciprocal lattice vectors have two important features: e*, e?,, 
■ ■ ■ , e* d are linearly independent; and fcj G Z, % — 1, 2, . . . , d. In the direct space, the lattice 
vectors have the same properties. 

The main idea of CAM is to use periodic structures to approximate quasicrystals. For 
a d-dimensional quasicrystals, its reciprocal lattice vectors k can be expressed by d linearly 
independent reciprocal vectors, e*, e^, • • • , e^, 

k = pie* + p 2 e* 2 + ■ ■ ■ + p d e* d . (9) 

It is important to note that the pi G TZ, that is, k cannot be represented by linear combi- 
nations of e* with integer-valued coefficients. However, the quasiperiodic function 0(r) can 
be expanded as in the following form 

0(r) = We HLkyr/L , r G [0, 2irL) d . (10) 

k 

If there exists a real number L such that Lpi G Z or Lpi can be made arbitrarily close to 
an integer for alH = 1, 2, . . . , d and k, which means 

\L-(p 1 ,p2,...,Pd)-([Lp 1 \,[Lp 2 \,...,[Lp d })\ l ~ ^0, for all k, (11) 
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where [•] rounds the number • to the nearest integer. Then numerical methods designed 
for periodic structures can be used to treat quasicrystals. Without loss of generality, we 
can always choose (1,0,..., 0) as one of the primitive reciprocal vectors. Therefore, the 
real number L becomes an integer. The problem of determining L is a well-known problem 
called the Simultaneous Diophantine Approximation (SDA) in number theory^. In view of 
numerical computability, the integer L must be as small as possible. 

From the above description, the condition of SDA must be satisfied for the implementa- 
tion of CAM. The integer L is dependent on these irrational numbers Pi due to rotational 
symmetry and the desired precision of the approximation. In the subsequent numerical 



examples (in Sec. II C), we will find that L increases very quickly as the desired precision be- 



comes higher. This greatly increases the computational cost. For CAM, the computational 



box should satisfy the condition (11). Therefore, the edge length of computational box 



is close to -Dk = L x 2tt, or D = n ■ L x 2tt, n G M + . 



2. An Example of 12-fold Rotational Symmetry 



(-cosO/3), sin(7i/3)) e 2 A (0 ' !) (cos(7i/3),sin(7i/3)) 



(-cos(ti/6), sin(7i/6)) 

(-1, 0) < 

(-cos(7t/3), -sin(7i/3)) 

(-cos(7i/6), -sin(7i/6)) 




(0, -1) 



(cos(7t/6),sin(7i/6)) 

> e !*=(b0) 

(cos(7r/6), -sin(7r/6)) 
(cos(7t/3), -sin(7i/3)) 



FIG. 2: Sets of reciprocal lattice vectors {k} of 12-fold rotational symmetry in 

2-dimensional space. 



In order to make the above method clearer, we take an example of 12-fold rotational 
symmetry in two-dimensions to demonstrate this point. A diagram in the reciprocal space 
is shown in Fig.[2] Let two noncollinear vectors = (1,0) and = (0,1) comprise the 
primitive reciprocal vectors, other reciprocal vectors of the 12-fold case are represented as 
linear combinations of and with real number coefficients, as shown in Fig.[2j If other 
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two noncollinear reciprocal vectors are chosen as primitive reciprocal vectors, similar results 
will emerge as well. The distinct nonzero coefficients for these 12 vectors are 1, 1/2, and 
v3/2. When applying CAM to this example of 12-fold vectors, we need an integer L such 
that 



L- 1 



V3 1 

~2~' 2 



L. 



v/3 



L 



2 J 



->• 0. 



(12) 



5. Applying CAM to Lifshitz-Petrich Method 



In CAM, we indeed use a crystalline approximant to approximate the quasicrystal. By 
substituting Eqn ( JlO) ) into Q, the Lifshitz-Petrich free energy density functional becomes 

F[<P\ = \ E b(l-|k| 2 ) 2 (g 2 -|k| 2 ) 2 -^(ki)0(k 2 ) 



ki+k 2 =0 

a 
~ 3 



E 0(ki)0(k 2 )0(k 3 ) + ^ E 0(ki)0(k 2 )0(k 3 )0(k 4 ). (13) 

k 1 +k 2 +k 3 =0 ki+k 2 +k 3 +k 4 =0 

Then we use semi-implicit scheme Q to minimize the free energy density functional. For 
CAM, the semi-implicit method becomes 



+ c(l-k 2 )V 



At 



k 2 ) 2 ] <W(k) 



At +e)0 4 (k) + a(0 2 )(k)-(0?)(k), (11) 



where (0™)(k) = J <ir0 m (r)e jkr , m 



2, 3. The right terms of the dynamics can be 
efficiently calculated by the pseudospectral method^. The Laplacian terms are computed 
in d- dimensional reciprocal space easily, while the convolutions can be calculated efficiently 
by Fast-Fourier transformation (FFT). Then the computational complexity is 0(iV log 
at each time step, iV is the number of degrees of freedom. 



B. Projection Method (PM) 

1. Method Description of PM 

As mentioned above, CAM computes the crystalline approximants. However, the approx- 
imation method cannot evaluate the free energy density exactly because of the approximate 
error of SDA. Therefore it is necessary to design a new numerical method that improves 
the calculation of quasicrystals by avoiding the approximate error of SDA. From Fourier- 
space crystallography^ 7 * 19 * 20 !, there exists an equivalent n-dimensional representation of a 
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<i-dimensional quasicrystal (n > d). In the higher-dimensional description* 17 * 19 * 2 ^ , the recip- 
rocal vector k of a (i-dimensional quasicrystal is, 

k = hipl + h 2 p* 2 H h h n p* n , hi G Z, (15) 

with vectors p* G TZ d of rank n. In the following, one should keep in mind that k G lZ d . We 
can redefine the selected reciprocal vectors p* with n components. Then the quasicrystal 
can be represented by these redefined basic reciprocal vectors with integral coefficients in 
n-dimensional space. This representation has the advantage, for our purpose, that the d- 
dimensional quasicrystal is periodic in n-dimensions. In order to make the description more 
clear, we introduce the n-dimensional reciprocal lattice. Assume that n-dimensional vectors, 
b*, t>2, . . . , b* , are the primitive reciprocal vectors of a 1st Brillouin zone in n-dimensional 
reciprocal space, the reciprocal vector of n-dimensional periodic structure can be expressed 
as 

H = h x h\ + h 2 b* 2 + ■■■ + h d b* n , (16) 

where the coefficient hi G Z, and H G lZ n . The correspondingly primitive vectors of the 
direct space hi satisfy the dual relationship (JBJ. We want to use the n-dimensional reciprocal 
vector H to represent the (i-dimensional quasicrystal. Then we can solve a (i-dimensional 
quasicrystal as a periodic structure in n-dimensional space. The key point of implementing 
the above idea is to provide an operator to project the n-dimensional structure into d- 
dimensional space. 

In order to solve the problem, we propose a novel Fourier expansion for the (i-dimensional 
quasiperiodic function acting <i-dimensions 

g(r) = Y J g(ny [is - n)T - r] . (17) 

H 

In the expression, r G TZ d , H G TZ n , and S is the projective matrix which connects the 
(i-dimensional physical space with the n-dimensional reciprocal space. We note that the two 



representations (15) and (16) can be used to describe the same quasicrystal. The recip- 
rocal vectors of a (i-dimensional quasicrystal can be represented by (i-dimensional recipro- 
cal vectors p* with integral coefficients, and also the integral combinations of extended n- 
dimensional primitive reciprocal vectors b*. Therefore, we can obtain the projective matrix 
S through projecting the n-dimensional reciprocal vectors bj into (i-dimensional reciprocal 



space, i.e., Pi = (S ■ b)», % = 1, 2, . . . , n. The j-th component of the projected reciprocal 
vector p* can be expressed by b* 

n 

Pij = s i mb *im> 3 = !> 2, • • • , d, (18) 

m=l 

where 6* m is the m-th component of the z-th n-dimensional reciprocal vector b*. Then these 
coefficients Sj m e 1Z, j = 1, 2, . . . , d, m — 1, 2, . . . , n, form the d x n-order nonzero projective 
matrix S. For the least computational expense, the dimension n of the extended space 
must be the smallest determined by the order of the elements of the symmetric group^^l. 
For example, 5-, 8-, 10-, and 12-fold symmetric quasicrystals, the dimension of embedded 
space is four. However, 7-, 9-, and 18-fold symmetric quasicrystals must be restricted to 
six-dimensional space. The projective matrix S is not unique, which is determined by the 
selection of reciprocal vectors b*. The represented coefficients of reciprocal vectors are 
dependent on the selection of primitive reciprocal vectors, however, the reciprocal vectors of 
a quasicrystal are unique. Therefore, the selection of basic reciprocal vectors b* as well as 
the projective matrix S is irrelevant to the quasicrystal. Considering a periodic structure, 
the projective matrix S degenerates toadx <i-order unit matrix. 

Furthermore, PM can be extended to calculate one-dimensional quasicrystals even though 
the notion of rotational symmetry does not exist in one dimension. If an energy functional 
has one-dimensional quasiperiodic structures with different incommensurate scales, for ex- 
ample, 

g{x) = C sin(x) + C\ sin(gix) + C 2 sm(q 2 x) H , (19) 

where the common multiples of 1, qi, q 2 , ■ ■ ■ , are irrational numbers in pairs, the projective 
matrix becomes a vector 

S = (l, gi ,q 2 ,...). (20) 

Then our proposed method can be applied to calculate one-dimensional quasicrystals. 
In the following, a lemma is given to indicate which variable should be computed in PM. 



Lemma 1. For a d- dimensional quasiperiodic function g(r), under the expansion p7y , we 
have 

lim i fdTg{T)=g{K) . (21) 

V^oc V J H=0 
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Proof. Firstly, we note that 

lim I- /drexp((<S-H)-r) = 5(<S-H). (22) 

V^oo V J <■ J 



Therefore, 



lim - / drg(r) = lim - / dr Vg(H)e l[(5 ' H)T - rl = o(H) . (23) 

j H 



Then we just need to prove that the expressions (21) and (23) are equivalent. From the 



definition of (16), n-dimensional reciprocal vector H is the integer- valued combinations of 
linearly independent primitive reciprocal vectors b*, i — 1, 2, . . . , n. Because the projective 
matrix S is linear, the d- dimensional projected vector S • H is also the integer-valued com- 
binations of the (i-dimensional projected vectors p* = S ■ b*. S ■ H = is equivalent to the 



integral coefficients hi in Eqn. (23) to be zero. It means that the n-dimensional reciprocal 



vector H = 0. □ 

Remark 1. From Lemma 1, the Fourier coefficients (?(H) rather than g(S ■ H) should be 
computed in PM. 

In PM, a quasicrystal is computed in n-dimensional reciprocal space as a periodic struc- 
ture, then the n-dimensional structure is projected into ci-dimensional space to obtain the 
(i-dimensional quasicrystal. Since different periodic phases have their own periodicity, the 
appropriate computational box is important to determine the final morphology of solutions, 
especially for complex phases. For example, in diblock copolymer systems^, the lamellae 
phase can be obtained easily in any computational region, however, for complex gyroid pat- 
tern, the computational box should be close to its period. For more complex quasicrystal, 
the computational box should be estimated carefully before computing. We also note that 
a equilibrium periodic structure is not only the minimum of a free energy density F with 
respect to order parameters, but also with respect to the unit cell^. Therefore, the unit 
cell in n-dimensional space, B = [bi,b 2 , . . . , b n ], should satisfy dF/dB = 0. We will give 
a estimation formula for the unit cell of quasicrystals in the Lifshitz-Petrich model (see 



Sec.IIB3) 



The PM is just with the help of n-dimensional reciprocal space to compute the d- 
dimensional quasicrystal. The physical space variable r always belongs to (i-dimensional 
space. Therefore an energy functional including (i-dimensional quasicrystals should not be 
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up to describing n-dimensional structures. The PM is able to compute quasicrystals rather 
than crystalline approximants. The PM is also using the periodic condition conveniently in 
higher- dimensional reciprocal space . Compared with CAM, the most important advantage 
of PM is that the method overcomes the restriction of SDA. Therefore PM can compute 
the free energy density to high accuracy numerically. Meanwhile, the proposed method can 
implement in a n-dimensional reciprocal unit cell to reduce computational cost. 



2. An Example of 12- fold Rotational Symmetry 



As an example, we take the 12-fold rotational symmetry in two dimensions to illustrate 
the PM. As Fig. [2] shows, the 12 reciprocal lattice vectors cannot be represented by two 



b 4 =(0, 0, 0, 1) 




b 3 = (0, 0, 1, 0) 



(-1,0,0,0) ^ 



b = (0, 1, 0, 0) 



> b;=(i,o,o,o) 



(0, 1, 0, -1) 
(1,0,-1, 0) 



(0, 0, 0, -1) 



FIG. 3: Sets of reciprocal lattice vectors {H} of the 12-fold rotational symmetry with 4 

components. 



noncollinear vectors and with integral coefficients. However, as Fig. [3] shows, if the 
vectors of b^, h^, bg, and b^ with 4 components are chosen as primitive reciprocal vectors, 
other reciprocal lattice vectors can be represented as integral combinations of these primitive 
vectors. The selected 4 vectors make up a basis in 4-dimensional space. The 2-dimensional 
12-fold example is related to a periodic structure in 4-dimensional reciprocal space. The 
projected vectors p* = «S-b*, 2 = 1,2,3,4 also can be expressed by e*, in 2-dimensional 
space, i.e., 

(i — lV „ „ . (i - 1)-7T „ , nA . 
Pn = cos e n , p i2 = sin e 22 . (24) 
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Therefore, the projective matrix is 

, 1 cos(7r/6) cos(7r/3) . 
S= ( V 1 ' V 1 ! | . (25) 

sin(7r/6) sin(7r/3) 1 

3. Applying PM to Lifshitz-Petrich Method 

Applying PM to the Lifshitz-Petrich model, the order parameter <p(x, y) is expanded as 
0(x, y) = <K H ) ex P ■ H ) T ' K vf}, (26) 

H 

where H is a n-dimensional vector, S is a 2 x n-order projective matrix. For quasicrystals, 
the system volume V should go to infinity. Based on Lemma 1, the Lifshitz-Petrich free 
energy density functional becomes 

FW>v)\ = \ E W 1 - + aDfh 2 - (a! + gl)f - e}#(H0 0(h 2 ) 

Hi+H 2 =0 



E 0(H!)0(H 2 )0(H 3 ) + ^ E 0(H X )0(H 2 )0(H3)0(H 



3 ^ rv /rv /rv u/ 4 

H!+H2+H 3 =0 Hi+H 2 +H 3 +H 4 =0 



4 



(27) 



where #i and g 2 are defined by 



«s * h = ( E sh E ^i 6 i*' E S2j E h 3 b s*) = fi,2 ) T ' ( 28 ) 

8=1 j = l j=l j = l 

/ij G 2, Sij and s 2 j are the components of the projective matrix S, and bji is the i-th 
component of the primitive reciprocal vector bj. The n-dimensional Fourier coefficient 0(H) 
can be solved by the semi-implicit method ^ 



■At 



At 

where the quadratic term and the third term are 



' e)0 4 (H) + a(0 2 )(H)-(^)( H ), (29) 



)(H)= ^(Hi)0t(H 2 ), ($)(H) = E &(Hi)&(H a )# t (H 8 ). (30) 

Hi+H2=H Hi+H2+H3=H 



In Eqn. (29), the linear term can be solved easily. The nonlinear terms are n-dimensional 



convolutions in reciprocal space. Directly computing will result in expansive computational 
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cost. However these convolutions are local point multiplication in n-dimensional direct 
space. Therefore we use the pseudospectral method^ to treat these terms by FFT. It 
should be emphasized that the PM is implemented in n-dimensional reciprocal space than 
(i-dimensional physical space. Computing these convolutions in n-dimensional direct space 
is to reduce computational complexity. 

For PM, we give a method to estimate the n-dimensional unit cell for the Lifshitz-Petrich 
model^. Without loss of generality, we can choose a proper coordinate system such that 
bu 7^ 0, bij = 0, when j > i. If B = [b!,b 2 ,. . . , b n ] is a n-dimensional unit cell of a d- 
dimensional quasicrystal. The first deviations of the free energy functional with respect to 
bij should be zero, where i — 1, 2, . . . n, j — 1, 2, . . . , i, i.e., 

J2 h ^9isu + - (gl + giW - {gl + gl)][i + g 2 - 2(0? + gl)} ■ l<KH)| 2 = o. (31) 

H 



Then the unit cell can be obtained by solving the Eqn. (31). Since the Fourier coefficients 
0h are unknown beforehand, in practice, only primary reciprocal vectors with equal Fourier 
coefficients are considered in estimating the unit cell. 



C. Computational Complexity 

In this section, we will give a general analysis on the computational complexity of CAM 
and PM in solving the two-dimensional Lifshitz-Petrich model. The computational complex- 



ity of CAM is dependent on the approximate accuracy of SDA ( 11 ) and the numerical resolu- 
tion. However, the PM overcomes the restriction of SDA, whose computational complexity 
is only dependent on the numerical resolution. As described in Sec |II A i] and Sec. II B 1 , both 



CAM and PM can use pseudospectral method in computing the two-dimensional Lifshitz- 
Petrich model. Then the number of basic functions used in the reciprocal space is equivalent 
to the number of discretized points in the direct space. In order to guarantee the same nu- 
merical resolution, we assume that the mesh step size is Ax both for CAM (2 dimensions) 
and PM (n dimensions) in the direct space. 

In the Lifshitz-Petrich model, the CAM is implemented in two dimensions. At a desired 
approximate error Esda as well as the integer L, the grid size is x N^, 



Ac 

Ax 



D k = 2tt x L, (32) 
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where [•] rounds the number • to the nearest integer. Then the computational complexity of 
CAM 



0(N tM -2N*logN k ) 



(33) 



where iV^k is the number of time iterations. 

In PM, the 2-dimensional quasicrystal is represented in n-dimensional reciprocal space 
as a periodic structure. The dimension of n is dependent on the rotational symmetry of a 
quasicrystal. Therefore the 2-dimensional quasicrystal can be computed in a n-dimensional 
unit cell. We assume that the size of the n-dimensional unit cell is Dh hi the direct space. 



Then the number of the basic functions is iVjj, with 



iVi 



H 



D 



H 



Ax 



(34) 



As discussed in Sec. II B 1 the pseudospectral method in PM is used to solve the dynam- 



ical equation (29) in n-dimensional reciprocal space with the help of FFT. Therefore the 



computational complexity of PM is 



0(N t , n -nN£\ogN H ), 



(35) 



where N tt n is the number of time iterations. 



From the estimation formulas (33), the computational complexity of CAM is a function 



of the integer L, which is dependent upon the approximate error of SDA. In next section, 
we will find that 3> Dh as the approximate error of SDA decreases. The computational 
complexity of CAM may be larger than that of PM in situation where higher numerical 
accuracy is required. More importantly, PM can compute quasicrystals and their free energy 
density to high accuracy without any approximate error of SDA. However, CAM always has 
the approximate error of SDA unless the computational box goes to infinity which results 
in unaccepted computational cost. 



III. NUMERICAL RESULTS AND DISCUSSION 

We will demonstrate the behavior of the two numerical methods, the CAM and the PM, 
based on the two-dimensional Lifshitz-Petrich model. In previous research^ 22 * 26 *, the Lifshitz- 
Petrich model showed that if q is chosen around 2 cos(7r/12) one can obtain a 2-dimensional 
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quasicrystal with dodecagonal (12- fold) symmetry, the earlier work found that no choice of 
q yields globally stable octagonal or decagonal symmetric pattern. However, in the recent 
study, the decagonal quasicrystal is found to be stable in the Lifshitz-Petrich model^. Here 
we only consider the dodecagonal symmetric structure. 



A. Computational Complexity of Dodecagonal Symmetric Structure 

For the dodecagonal quasicrystal (DDQC), the initial reciprocal vectors at which the 
Fourier coefficients are nonzero, for DDQC in the Lifshitz-Petrich model is shown in Fig. [4], 
which contain two 12-fold stars of wave vectors, one on |k| = 1 and other on |k| = q = 
2 cos(7r/12)^. The specific reciprocal vectors are given in Tab.| when e* = (1, 0) and 



ML 
i 



\ A 

\ i 
\ , t 



FIG. 4: Initial reciprocal lattice vectors for DDQC. 



TABLE I: Initial reciprocal lattice vectors for DDQC in 2-dimensional space with 

q = 2cos(vr/12). 



(cos(jvr/ 6 ) ! sin (W 6 ))> J = 0, 1, . . . , 11 
(q cos(j7r/6 + vr/12), q sin(jV/6 + vr/12)), j = 0, 1, . . . , 11 



e* 2 = (0,1), as shown in Fig[2j are chosen as the basic reciprocal vectors in 2 dimensions. 
The distinct nonzero represented coefficients of these initial reciprocal vectors for DDQC 
phase are 1, 1/2, a/3/2, 2cos(7r/12), cos(7r/12) and \/3 cos(7r/12). The condition of SDA 



( 11 ) must be satisfied when one uses CAM. Fig. 6 shows the trend of the approximate error of 
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x 10 

FIG. 5: The trend of the approximate error of SDA, Esda-, as a function of the integer L. 

SDA, Esda, as the integer L increases. From the Fig.[6j one can find that the convergence 
of Esda is not uniform. Tab. [IT] gives the minimal integer for desired approximate error 



TABLE II: The minimal integer L for desired approximate error of SDA, Esda- 



Esda 


0.19098 


0.17486 


0.07042 


0.04953 


0.03583 


0.02961 


0.01936 


L 


30 


208 


410 


3404 


6016 


32312 


82262 



The size of computational box used in CAM is = 2ttL. 



Esda, and the corresponding size of computational box used in CAM. Then the initial 
reciprocal vectors in CAM are [L ■ k]. From the table, we find that L increases quickly 
as the approximate error becomes small. Accordingly, the computational cost will increase 
greatly. We also find the Esda ~ 0.19098 (L = 30) is the least requirement for computing 
dodecagonal symmetric structure, it is consistent with the result in Ref. [22]. The initial 
reciprocal vectors in 4-dimensional space representation of PM are shown in Tab. |III| when 
the primitive reciprocal vectors are b{ = (1,0,0,0), = (0,1,0,0), bg = (0,0,1,0), and 
b 



to 



, A — (0, 0, 0, 1), as shown in Fig|3j We can use the 24 initial reciprocal vectors in Tab. Ill 
estimate the unit cell of DDQC in 4-dimensional space. The projective matrix S of PM is 
given by Eqn. (25). From the estimation formula (31) and the dual relationship ([6]), the size 
of the unit cell in n-dimensional direct space is Dh = 2-k. 

In the following, we will compare the computational complexity between CAM of different 



unit cell with D k and PM as discussed in Sec. |II C| In each time step, for the dodecagonal 
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TABLE III: Initial reciprocal lattice vectors {H} for DDQC in PM with q = 2cos(7r/12). 



|5-H| = 1 


(0 1 -1) (0 -1 1) (1 0) (-1 0) (0 1 0) (0 -1 0) 
(0 1 0) (0 -1 0) (0 1) (0 -1) (-1 1 0) (1 -1 0) 


|5-H| = q 


(1 1 -1) (-1 -1 1) (1 1 0) (-1 -1 0) (0 1 1 0) (0 -1 -1 0) 
(0 1 1) (0 -1 -1) (-1 1 1) (1 -1 -1) (-1 -1 1 1) (1 1 -1 -1) 


The projective matrix S is given by Eqn. (25 


). 



symmetric structure, the computational complexity of CAM is 



(36) 



and the computational complexity of PM is 



C H = O (4A£l gjV H ) =0(4 (J?-} lo 




(37) 



For L = 30 with Esda ~ 0.19098 in CAM, the computational complexity Ck < Ch when Ax 
is smaller than 0.20943951. However, for L = 208, C k < C H only if Ax > 500.6549, which 
is much larger than the size of the unit cell Dn. It cannot be implemented numerically. 
The computational complexity Ch of PM is always smaller than that of Ck- For higher 
approximate accuracy of SDA with larger L, Ch < Ck, and total computational complexity 
of PM may be less than CAM. 

Then we compare the total computational complexity including time iterations of CAM 
(with L = 30) and PM through numerical experiments. All the methods considered here 
have been implemented in C language. Fourier transforms are computed using the FFTW^l 
package. The codes were run in the same dedicated computer, a Intel(R) Xeon(R) CPU 
E5450 @3.00GH memory 16 G under linux. The time step size, At, is always selected as 0.1 
for both methods. To measure error we use l°° norm, 



fcAM = max 

k 



(k) 



(38) 



for CAM, and 



f PM = max { j (H) 
for PM. In the following simulations, the error tolerence is set as 1.0 x 10~ 8 . 
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(39) 
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F CAM = -3-457837609e-03 



F D . . = -3.524067379e-03 

PM 




0.2094 0.1795 0.1571 0.15 



FIG. 6: The total CPU time, the free energy density and the number of time step at the 
error 10~ 8 of CAM (L = 30) and PM, with different mesh step size Ax at a set of 

parameters: c = 50, e = 0.015, a = 1.0. 



Fig. [6] shows the total CPU time, the free energy density, the number of time step at the 
desired error 10" 8 , of CAM (L = 30) and PM at a set of parameters, c = 50, e = 0.015, 
a = 1.0, with different discrete resolution. The number of time step of CAM on different 
grids is same, N t ^ = 848. PM also has the same time step to achieve the the error 10~ 8 on 
different grids, N t ^ = 472. It may be that the accuracy of spacial discretization is enough 
for the calculated structure at the set of parameter. For the same reason, the free energy 
density calculated by PM is nearly the same on these grids, while the CAM has the equal 
free energy density. However, the number of time step of CAM, N t ^ is alway large than that 



of PM, iV t) H- As discussed in Sec. Ill B it may be that PM keeps the rotational symmetry, 
while CAM dose not. Then PM can reach the same accuracy quickly. The free energy 
density computed by PM is lower than the CAM. It may result from the approximate error 
of SDA in CAM. As Fig[6] shown, the cost of CPU time of PM is lower than that of CAM 
until the mesh step size Ax is less than about 0.15. It is different from the above analysis 
about the computational complexity in each time step because of the different time step in 
PM and CAM. 

From these results, we find that with relative large mesh step size Ax, PM can obtain 
enough numerical accuracy with less computational cost than CAM (L = 30). CAM always 
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has approximant error unless — > oo. And as discussed in the section, higher accuracy 
of approximant error will result in heavy computational burder. Therefore, PM has less 
computational complexity than CAM to high accuracy in computing these quasicrystals 
represented in 4 dimensions, such as dodecagonal symmetric structures. 



B. Dodecagonal Symmetric Structure 



The initial reciprocal coefficients, as mentioned in Sec. HI A , at which the Fourier co- 
efficients are nonzero are used as initial values to find equilibrium dodecagonal symmetric 
structures. Unless otherwise specified, the simulations of CAM are performed on a 256 x 256 
grids with L = 30, and the simulations of PM are performed with iVjj = 32 4 plane waves in 



the following. The projective matrix S of PM is given by Eqn. (25). We also find that the 



morphologies and the free energy density will not change with denser grids. The morphology 
of the dodecagonal approximant computed by CAM in physical space is given in Fig. 7(a)[ a) . 
It is a periodic structure. The DDQC calculated by PM is shown in Fig. 7(b)[ b), which can 
be viewed as a part of the dodecagonal approximant. Both morphologies demonstrate the 





(a)Dodecagonal crystalline approximant (b)DDQC 

FIG. 7: The morphologies of (a). Dodecagonal approximant computed by CAM and (b). 
DDQC calculated by PM at e = 0.015, a = 1, c = 50, q = 2 cos(tt/12). 

12-fold orientational symmetry in the physical space, at least locally. However, further anal- 
ysis about the relationship of the Fourier coefficients will come to a different conclusion. 
From the rotational symmetry, each 12 reciprocal vectors with dodecagonal rotational sym- 
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metry on the ring |k| = 1 (|«S • H| = 1), or |k| = q (\S ■ H| = q), should have equal Fourier 



coefficients. Tab. IV shows the principle Fourier coefficients calculated by CAM. Tab.fV] 
gives these Fourier coefficients on principle reciprocal vectors calculated by PM. From these 
results, we can find the 12 Fourier coefficients calculated by CAM on each characteristic ring 
are not equal. In contrast, as Tab. [V] shown, PM can obtain the 12 equal Fourier coefficients 

TABLE IV: The Fourier coefficients of the principle reciprocal vectors for dodecagonal 
crystalline approximant computed by CAM with L = 30 at e = 0.015, a = 1, c = 150, 

q = 2cos(vr/12). 



|k| = L 


(30,0) (-30,0) (0,30) (0,-30) 


6.106618210e-02 


(26,15) (-26,-15) (15,26) (-15,-26) 
(-15,26) (15,-26) (-26,15) (26,-15) 


5.928204525e-02 


\k\ = [L-q] 


(41,41) (-41,-41) (-41,41) (41,-41) 


5.458000478e-02 


(56,15) (-56,-15) (15,56) (-15,-56) 
(-15,56) (15,-56) (-56,15) (56,-15) 


5.683018766e-02 



on each ring. We also find the same phenomenon at different parameter coordinates. For 

TABLE V: The Fourier coefficients of the principle reciprocal vectors for DDQC 
computed PM at e = 0.015, a = 1, c = 150, q = 2cos(vr/12). 



|5-k| = 1 


(0 1 -1) (0 -1 1) (1 0) (-1 0) 
(0 1 0) (0 -1 0) (0 1 0) (0 -1 0) 
(0 1) (0 -1) (-1 1 0) (1 -1 0) 


5.856822141e-02 


\S-k\=q 


(1 1 -1) (-1 -1 1) (1 1 0) (-1 -1 0) 
(0 1 1 0) (0 -1 -1 0) (0 1 1) (0 -1 -1) 
(-1 1 1) (1 -1 -1) (-1 -1 1 1) (1 1 -1 -1) 


5.855442187e-02 



CAM, we also use higher approximate accuracy E$da ~ 0.07042, with the computational 
box Dk x Dk, = 8207T. The numerical experiments are implemented on a 4096 x 4096 
grids. We find that CAM can capture the dodecagonal approximant, however, it cannot 
obtain 12 equal Fourier coefficients on each ring as well as the above numerical experiments. 
Therefore, we come to a conclusion that PM can keep the non-crystallographic symmetry 



21 



to high accuracy, while CAM cannot. 



C. Free Energy Density 

We compare the free energy computed by CAM, PM and single-wave approximation 
(SWA). The SWA uses these principal Fourier vectors to approximate the free energy func- 
tional 29 . It translates the energy functional into a single- variable or mult i- variable function. 
Then the minimum of free energy can be obtained by minimizing the reduced energy function 
with these variables. The principal Fourier vectors of DDQC have been shown in Fig|4} The 
reduced free energy function of DDQC computed by SWA in the Lifshitz-Petrich model^EH 
is (C -> oo) 



12 



99 ( 



b{ + 0j) + 144(0? + 0?0 9 ) + 36O0 2 2 - 24(0 2 0, + 2 g ) - 8(0? + - 6e(0 2 + 2 ,)/a 2 



(40) 



0i,0g G TZ stand for Fourier coefficients on the |k| = 1 and |k| = q rings, respectively. The 
approximated minimum of DDQC can be obtained by minimizing the above equation with 
respect to 0i and <p q . The result of SWA can be viewed as a approximate solution of DDQC 



,x 10 




x10" 



(a) Free energy density 



200 




200 



(b)Interfacial energy (Laplacian terms) 



FIG. 8: (a). Free energy density calculated by PM and CAM, as a function of the penalty 
factor c, relative to that of SWA (c — >■ oo), with q = 2cos(7r/12), a = 1.0, e = 0.015. (b). 
Corresponding interfacial energy computed by PM and CAM. 



when c — > oo. Fig. 8(a) shows the free energy density calculated by PM and CAM, as a 
function of the penalty factor c, relative to that of SWA (c — > oo) with q = 2cos(7r/12), 
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a = 1.0, £ = 0.015. The free energy density Fqam computed by CAM is less than Fswa 
obtained by SWA until the penalty factor c increases to about 10. Then F C am is larger than 
Fswa as c increases. The reason is that the CAM cannot control the principal reciprocal 
vectors located on |k| = 1 and |k| = q when the penalty factor c increases, as shown in 



Fig. 8(b) In other words, CAM cannot calculate the interfacial energy (Laplacian terms) 
accurately when c is large. However the role of the differential terms in the Lifshitz-Petrich 
model is to keep the interactions at two characteristic scales. It is one of the reasons that 
the Lifshitz-Petrich model is able to stabilize quasicrystals^l. Correspondingly, the F C am 
diverges from the true free energy density. We have also used L = 410 which improves the 
approximate accuracy of SDA in CAM to observe the free energy density. It costed much 
computational amount on a 4096 x 4096 grids, however, the same problem appeared. In 
contrast, the free energy density Fpm calculated by PM is always smaller than Fswa for all 
c and converges to the Fswa as c — > oo. It is because the more basic functions are used in 
simulations by PM than that of SWA. And the approximated space of PM is more precise 
than that of SWA. Meanwhile, PM can maintain its principle reciprocal vectors located on 
scales 1 and q which is consistent with the model, as shown in Fig. 8(b)| Therefore PM can 
compute the free energy density for all parameters of high numerical accuracy. 



IV. CONCLUSIONS AND OUTLOOK 



In the article, we summary the CAM, and point out the advantage of using the periodic 
condition and the restriction condition of the SDA in this method. Then we propose a new 
numerical method, PM, which is implemented in higher-dimensional reciprocal space. The 
developed method overcomes the restriction of SDA and enables us view a quasicrystal as a 
periodic structure, and uses the periodic condition conveniently. The projection method can 
reduce the computational effort efficiently by computed quasicrystals in a higher-dimensional 
unit cell and using the pseudospectral method. By applying two methods to the Lifshitz- 
Petrich model, we analyze the computational complexity. The computational complexity 
of CAM is dependent on the approximate error and the numerical resolution, while that of 
PM is only dependent upon the numerical resolution. Specially, in computing dodecagonal 
symmetric pattern, PM has less computational complexity than that of CAM to high accu- 
rate. We also find that our approach can keep the rotational symmetry accurately, and more 
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significantly, the algorithm can calculate the free energy density without any approximate 
error of SDA. However, the dimensions of our computational space are dependent on the 
rotational symmetry of quasicrystals. For quasiperiodic structures with 5-, 8-, 10- and 12- 
fold symmetry, the computational space dimension is four. PM has accepted computational 
complexity. For a quasicrystal with of 7-, 9-, and 18-fold quasicrystals, the dimension of 
computational space is up to six. Therefore, our future work will focus on the improvement 
of the computational efficiency of our proposed method. Applying the proposed method to 
study real physical systems is also one of the major tasks in our further research. Finally, 
we should point out that PM can be also applied to study general <i-dimensional aperiodic 
structures^ whose Fourier spectrum consists of 5-peaks, H = Y^=i h%b*, hi G Z, of rank n 
(n > d) with basis vectors b*, i = 1, 2, . . . , n. 
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